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Abstract 

We have studied the elementary excitations in a two component Bose-Einstein 

condensate. We concentrate on the breathing modes and find the elementary 

excitations to possess avoided crossings and regions of coalescing oscillations 

where both components of the condensates oscillate with the same frequency. 

For large repulsive interactions between the condensates, their oscillational 

modes tend to decouple due to decreased overlap. A thorough investigation 

of the eigenmodes near the avoided crossings are presented. 
03.75.Fi,05.30.Jp 
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I. INTRODUCTION 



Bose-Einstein condensation in trapped dilute gases differs from the bulk case in the 
lack of translational invariance [l|-|3|]. At zero temperature the ensuing spatially varying 
condensates can be calculated reliably from the Gross-Pitaevskii equations [|J. At finite 
temperatures the region close to the transition point has been treated by a Hartree-Fock 
scheme MM- 

In order to obtain the spectrum of the condensed gases a variety of methods have been 
applied, most of them utilizing a Bogoliubov-de Gennes diagonalization of the dynamics 
linearized around the ground state and the eigenmodes and their corresponding 

frequencies are well known ||T2HT6||. The calculations of the temperature dependence of the 
excitations have still not reached decisive conclusions [17|-pT|. 

A new physical situation emerged when it was found that two hyperfine components, 
\F = 2, m = 2 > and \F = l,m = — 1>, of Rb can be induced to condense concomittantly 



in a single trap [p2j] . The two condensates see different potentials, and they tend to appear 
displaced from each other thanks to the effect of gravity. This is however, a technical feature 
only and it is possible to make the condensates in a totally symmetric environment |SJ . In 
such a situation, we have shown |24| by a Hartree-Fock treatment that the symmetrically 
formed condensate may find a state of lower energy by breaking the symmetry. If the 
experimentally realized condensate falls into the parameter range where this ooccurs remains 
to be seen. 

In this paper we discuss the oscillational spectrum of a symmetric double condensate. 
Other works investigating the low lying excitations of the double condensate use a variational 



approach |25[ and a direct numerical treatment p6|-|28[|. We linearize the Gross-Pitaevskii 
equations of the double condensate and solve for the eigenmodes. Because of the numerical 
burden of the computations, we only look at the radially symmetric (7 = 0) oscillations, the 
breathing modes JT2]JT^] . We have used time dependent integration methods to check the 
results of our eigenmode analysis, the details will be explained below. As results we have 
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found the oscillation modes to display avoided crossings as expected of quantum energy 
levels. 

The organization of the present paper is as follows. Section II investigates the small 
oscillations of the coupled Gross-Pitaevskii equations for the double condensate. Their dy- 
namics is linearized and the usual Bogoliubov-de Gennes expansion |H],§2|] is generalized to 
the two-component condensate case, which gives an eigenvalue problem, which is nonher- 
mitian, but has real eigenvalues because of the time reversal symmetry of the problem. In 
Sec. Ill we explain how to solve this numerically, and present methods to verify the results 
by direct time integration of the equations when perturbed from equilibrium. The results of 
the calculations are presented and discussed in Sec. IV. The main conlusions are presented 
in Sec. V. 



II. FORMULATING THE PROBLEM 

We consider a two-component Bose condensed gas in external harmonic potentials 

V l (r) = ^mny 2=1,2 (1) 

with the trap frequencies Qi and assume a pair potential of the atom-atom interaction of 
the form Vi = Vi5(r — r') with 



A-nh 2 a; 



m 



1,2,3 (2) 



where a\ and a 2 stand for the inter species scattering lengths and 03 is the scattering length 
between the two different atoms, and m is the mass of the atoms which is taken to be the 
same for both species. The grand canonical Hamiltonian can then be written in the form 

r h 2 h 2 

H = J dr{*J(r)[-— V 2 + Vi(r) - Mi]*i(r) + $|(r)[-— V 2 + V 2 (r) - H^(r) 
+-«itf}(r)*I(r)tf 1 (r)tfi(r) + ^«a*S(r)^(r)* a (r)* 2 (r) 

+^M(r)^(r)^i(r)^ 2 (r)} (3) 
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where ^(r) is the field operator of the bosons and /i« is the chemical potential. To include 
this in the Hamiltonian only shifts the energies, and introduces the dependence on the 
particle number directly into the Hamiltonian. 

Using the canonical bosonic commutation relations between the fields {^i} we obtain 
the equations 

*ft-^(r) = h^V 2 + V x {t) + v 1 i>\{v)i> 1 {r) + v 3 ¥ 2 (r)$ 2 (r) - ^^(r) (4) 
F) h 2 

ih-i> 2 (v) = h^V 2 + V 2 (r) + v 2 ¥ 2 (r)if 2 (r) + i*#J(r)*i(r) - /x 2 ]* 2 (r). (5) 

Taking the expectation values and assuming the products to factorize, we obtain the 
time dependent two-component Gross-Pitaevskii equations 
d h 2 

i/i— *i(r) = h^V 2 + V x {r) + ^i(r)| 2 + ^ 2 (r)| 2 - ^i]*i(r) (6a) 
F) h 2 

ih-* 2 (r) = h^V 2 + V 2 (r) + ^ 2 |^ 2 (r)| 2 + ^i(r)| 2 - // 2 ]* 2 (r). (6b) 
where the chemical potentials are chosen such that 

J dr\*i(r)\ 2 = z = l,2. (7) 



When we set the time derivatives in Eqs. (paj) and (pb|) equal to zero, we obtain the 
Gross-Pitaevskii ground state functions {^f, ^2} which we assume to be chosen real. They 
solve the equations 

A*i*f (r) = [-t^V 2 + Vx(r) + ^?(r) 2 + v^ 9 2 (r) 2 ]*f (r) (8a) 
2m 

^ 9 2 (v) = [-|-V 2 + V 2 (r) + v 2 n(r) 2 + ^?(r) 2 ]M/|(r) . (8b) 
2m 

As is customary, we now expand the deviation from the ground state in the difference 

<^(r) = * f (r) - ttf(r) (i = 1, 2) (9) 

which still preserves the operator character of the problem. 

We can now proceed in two ways: One way substitutes the expansion @ into the Hamil- 
tonian (U) and retains the terms to second order only. The linear terms in are found 



to cancel because solves the stationary equations fl8a|) and (|3b]). The eigenmodes of the 
remaining Bogoliubov problem are then obtained by diagonalizing the quadratic part to 
obtain the diagonal Hamiltonian in term of creation and annihilation operators {aj,, a u } and 
{bl, b u }, which reduce to the eigenmodes of the separated condensates when the mode cou- 
pling v 3 vanishes. In terms of these operators, the Hamiltonian takes the form of decoupled 
linear oscillator relations. 

Alternatively, we may proceed from the equations of motion (P) and linearize these. 
Finding the eigenmodes, then consists in looking for the eigenfrequencies of this linear time 
dependent problem. In the present paper we choose to approach the problem in this way. 

From the two-component time dependent Gross-Pitaevskii equations (||) we get the lin- 
earized equations 

ih—svx = [h { h ] + 2v 1 vf + v^fs^i + « 3 (*f <j*x + ^m5% + m 9 2 5^ 2 ) (10) 

L/C 

= \H { ho + 2v 2 ^f}5^ 2 + v 2 ^ff6% + u 3 (*f + ^ 9 2 5^l + *f*f^x) (11) 

where 

flffo = -^V a + V < (r)-/i i . (12) 

In writing these equations, we have utilized the fact that the ground state functions \&f 
have been chosen real. The time dependent linear equations ( TO ) and (|TTJ) can be solved by 
generalizing the Bogoliubov transformation [|3^] to the two-condensate inhomogeneous case 

<Wi(r) = ^[a«(r)a,e- iE "^ + /3«(r)*aj; e iE ^] (13) 
S9 2 (r) = ^[a( 2 )(r)^ e -^*/ n + ^(r)^ 4 ^*/*], (14) 

V 

where the sum over ^ is a sum over elementary excitation modes of the system. Inserting this 



ansatz into (|I0~D and (111) we find the double condensate Bogoliubov-de Gennes equations 



H^o + 2v^f + v-M 2 ] o& + v 3 ^ 9 2 (a^ + (3^) + v^ff3^ = Ev0 £) ( 15a ) 
+ 2^!M/f + v^f 1 /?« + «s*f + /^ 2) ) + «x*f a« = (15b) 



hPo + 2t**f + «S*f] «1 2) + T*¥f*S(a?> + $>) + tto¥?0?> = E v o^ (15c) 

hSI + 2v 2 n 2 + v 3 n 2 } Pi 2) + v,nn{^ ] + $>) + v 2 *?aw = -e v $\ (isd) 

The ensuing eigenvalue problem is nonhermitian, but it has a symmetry equivalent with that 
in the case of a single condensate: Taking the complex conjugate, exchanging the a:s and 
/3:s and reversing the sign of the eigenvalue E u , we generate equations identical with ([Top. 



This time reversal operation guarantees that the eigenvalues are real and come in positive 
and negative pairs. There is also one pair of eigenvalues which are zero. These eigenvalues 
correspond to the Gross-Pitaevskii solutions, and are not in this case interesting since they 
do not represent any oscillations. They correspond to the zero frequency collective modes 
discussed for the single condensate by Lewenstein and You |33[ . 



III. NUMERICAL METHODS 



A. Solving the secular equations 



We now consider the case where we have spherically symmetric trapping potentials with 
Oi/0 2 = y/2. The coupled Eqs. (||) and fllED are then considerably simplified. In that 
case, the obvious and simplest way to solve the Bogoliubov-de Gennes equations ( |T5|) is to 
discretize the solutions and the derivatives. This method works well if the problem can be 
reduced to one dimension. Thus the radial part is solved for, which means that we can only 
see the so called breathing modes Hl2| , |l4| of the condensates. When discretizing the solutions 
(3^ and the derivatives we get an eigenvalue problem 



Mi 



-Mi 



\ 



M 2 V 2 ^2 

-v 3 ^m -v 3 *f*f -v 2 mf -M 2 



v 3 ^{m 9 2 v 3 m{m 9 2 M 2 v 2 mf 



E 



a® 



(16) 



with 
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Mi = H { S + 2v^f + v^f i ± j. (17) 

The eigenvalue problem then concerns a matrix of size (4n) x (4n), where n is the number 
of gridpoints used to build up the eigensolutions. 

The Bogoliubov-de Gennes equations can also be solved using a basis set method, where 
the solutions are expanded in some suitable orthonormal set of functions. In the case of 
harmonic traps, it is favourable to use the eigenstates of the harmonic oscillator. This gives 
an eigenvalue problem for the expansion coefficients instead of the more direct method of 
discretizing the solutions. Expanding the solutions in some basis set works particularly well 
when the condensates are separated by gravity. However, in the effective one- dimensional 
problem we consider, we found the grid method to be superior. In order to evaluate the 
matrix in Eq. ( |i6|) we need the solutions of the Gross-Pitaevskii equations (|8|). These were 
obtained with the method of steepest descent, which has been successfully used in earlier 
work P^ipDJ on the nonlinear Schrodinger equation. 



B. Time dependent methods 

The numerical diagonalization of the secular equation can be checked by a simple time 
dependent method. Choosing an initial state which is not an eigenstate of Eq. (|6a|) and (|6b|) 



is expected to produce oscillations in the densities. Following the time evolution of these 
at an arbitrary point in the condensate and fourier transforming the corresponding signal 
gives us the spectrum of the eigenoscillations. This crude fourier method, where an initial 
state which slightly differs from the ground state is chosen, works in the spirit of kicking 
the condensates and letting them ring. It gives us the spectrum in one go but it is not 
very accurate since we have used only 50 gridpoints in the densities. It is not obvious how 
to get a good spectrum with this method whithout greatly increasing the computing work. 
Increasing the grid to one hundred points has been found to give better agreement with the 
method in Sec. III. A, but the drawback is that the computing time becomes forty times 
longer. This method is thus found to be of limited accuracy and has to be considered as a 



check on the diagonalization procedure from Eq. (fig). 

Another approach to calculating the eigenvalue spectrum is to follow the response of the 



system to an external driving force This is most easily done directly on the nonlinear 
equations (|6|) treated as coupled scalar equations. As discussed in Ref. j3lfl , we add to the 
external trapping potentials Vi{r) the driving 

V d {r) = dcos(kr + u d t), (18) 

where uo^ is a variable driving frequency, and A; is a suitably chosen scaling parameter, whose 
value does not greatly affect the results. We can then follow the response of the condensate 
densities as functions of the driving frequency uid by fourier transforming the corresponding 
changes in the densities as explained above. This method is more accurate and also more 
transparent, since the condensate response is truly dramatic when Ud approaches one of the 
resonance frequencies. One problem is that with a sufficiently strong driving d, harmonic 
generation is also seen at multiples of the true resonance frequencies. The nonlinearities of 
the dynamic evolution mixes the frequencies. 



IV. RESULTS 

The results of our numerical calculations are reported in Fig. |l|. The continuous lines 
are the modes obtained from the secular equation ([16]) as explained above. Because of the 
spherical symmetry, the equations are reduced to a one- dimensional eigenvalue problem, 
which means that we see only the so called breathing modes of the condensates. At v 3 = 
the condensates do not interact with each other and behave like uncoupled systems. On the 
left hand side we number the eigenfrequencies as Jjl to fjlO starting from below. At v% — 
the frequencies fj2, f)4, Jj7 and jj9 belong to condensate one, the intertwined ones belong to 
condensate two. The spectrum is seen to display avoided crossings between Jjl and ft 2 or ft 6 
and ft 7 around v% = 0.07. 

The discrete points given in Fig.l derive from the response of the kicked condensate as 
described in Sec. III.B. As explained there, the results are not very accurate but serve as 
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an overall check on the spectrum. We found, that choosing a typical initial state to be the 
Gross- P it aevskii solution with t> 3 slightly shifted (At> 3 = 0.001) gives the best spectrum for 
levels up to Jj4. With this method it is possible to detect the four lowest frequencies without 
any major inaccuracy; for higher frequencies the agreement is less satisfactory. The solution 
of the time dependent Gross-Pitaevskii equations have to be computed for long times in 
order to get a high resolution spectrum. This was done with the grid method as explained 
earlier. The results from the fourier method shown in Figs.|||-f| used fifty gridpoints in the 
density. This gives stable iterations and allows us to solve for long times. The frequencies 
of the spectrum are slightly higher than the spectrum calculated from Eqs. (0). Fig. |5| 
shows a portion of the spectrum in Fig. [I] when we have increased the number of grid points 
to 100. We see that the frequencies obtained uniformly tend to approach the results from 
the secular equations. The fortyfold increase in computer time, however, makes it costly to 
cover the full spectral range shown in Fig. |l|. 

The eigenoscillations of the condensates are most clearly observed when the spectrum is 
investigated using the method of an external driving field. This method produces improved 
results at least for the three lowest frequencies. In this paper we have used d = 0.01 in 
order to keep Vd « Vi\^/f | 2 and set k = 1/2; see Ref. |JT|. In Fig. ^ we show the observed 



response of the two condensates at and near the lowest level with t> 3 = 0.07. To the left, we 
see the dramatic increase of the magnitude of the response when uj^ is changed from 3.4 to 
the resonance value 3.8. As this calculation has been done in the region near the avoided 
crossing, we find the same frequency nearly equally strong in the fourier transforms of either 
condensate (to the right in Fig. ||). In Fig. [7|, we show the radial densities of condensate 
one and two. We see that, already for t> 3 > 0.03, the repulsive interaction tends to separate 
the condensates thus decreasing their coupling. 

In order to illuminate the behaviour of the coupled double condensates, we consider the 
levels Jjl and Jj2 near t> 3 = 0.07, where an avoided crossing is suggested. Near t> 3 = 0.0, level 
ft 2 is an excitation in condensate 1 and fjl is condensate 2. If we look at the amplitudes 
ai,/?i,«2 and fa at t> 3 = 0.04, Fig. ^0 shows that the dominating oscillation amplitude 
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is still «2; the amplitudes ai and (3\ in Fig. §a are clearly smaller. If we compare these 
amplitudes with the densities of the condensates, Fig. [7|, we find that condensate 1 has no 
oscillational amplitudes after r « 2.5, because the density goes to zero. The condensate 2 
reaches out to about r ^ 3.5 because it is repelled by condensate 1. The features in the 
oscillation amplitudes ot\ and (3\ around r « 2 are caused by the coupling to condensate 2; 
for «2,/?2 the analogous effects are seen near r « 0. Beacuse the coefficients in Fig. || are 
eigenmodes of the coupled problem, there are areas where the components nearly decouple; 
thus near r m 2.2, the oscillations are nearly purely taking place in condensate 1. 

When we follow the oscillational level Jjl to V3 = 0.08, Fig. |9|, we assume that an 
avoided crossing has been passed. Thus the oscillational amplitudes have been transferred 
to condensate 1, which is verified in Fig. ^a. The amplitudes in condensate 2, Fig. ^|b, are 
clearly smaller. They disappear rapidly near r « 0, beacuse for this interaction strength, 
condensate 2 is pushed well away from the center of the trap, see Fig. [7|. Beacuse of the 
decreased overlap between the condensates, they influence each others oscillation amplitudes 
far less than for t> 3 = 0.04. 

Table I summarizes the results of all our numerical calculations. The frequencies of the 
five lowest eigenmodes are reported at V3 = 0.002 and V3 = 0.02. The first column is the 



result from the Bogoliubov-de Gennes secular equations fll5D . The next column labelled 
Fourieri is obtained from the spectral response of the kicked condensates and the column 
Fourier2 from the driven condensate response. Comparing the columns we obtain a good 
picture of the accuracy of our numerical methods. At V3 = 0.002 the frequencies agree well. 
At w 3 = 0.02 the frequencies calculated with the Fourier methods are less satisfactory as 
compared with the BdG-spectrum. Unfortunately, it is very easy to excite harmonic gener- 
ations in the vicinity of level jj4 and jj5, which makes it difficult to find the eigenfrequencies 
for levels higher than jj4. This may also explain the lack of improvement in Fourier 2 over 
Fourieri at these frequencies. 
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V. CONCLUSIONS 



We have generalized the Bogoliubov-de Gennes method to the case of two coupled con- 
densates. The ensuing eigenvalues are solved for the spherically symmetric excitations, the 
breathing modes. The calculations are numerically demanding, and obtaining the full spec- 
trum in this way seems to be beyond the numerical capacity of our approach. We have also 
started from a symmetric ground state, even if we know that the condensates may separate 
due to symmetry breaking This we have done mainly for numerical reasons, but our 
investigations of the condensate stability indicate that the range of interaction we are mainly 
interested in here, V3 ^ 0.03, may lead to stable condensates located on top of each other. 



We report on these results in a separate communication [[34 



We have used time dependent integrations to check the behaviour of the energy eigenval- 
ues in an indenpendent way. Starting from a perturbed state, the system evolution should 
contain all the eigenfrequencies of the system; here they are restricted to the symmetric ones 
of course. The advantage of this method is that a fourier transform brings out all frequen- 
cies at the same time. This method is found not to be very accurate. However, the method 
seems to verify the overall behaviour of the modes but in most regions the results fall above 
those based on the linearization method. This is, however, not true in all parameter ranges, 
see Fig. [TJ. A more exact, but quite time consuming method brings the results of the time 
integration closer to those derived from the linearized method (see Fig||]). 

By perturbing the condensates at a single frequency, the resulting response grows dra- 
matically at the resonances (see Fig. |6|), but to obtain the full spectrum the calculations 
have to be repeated for all frequencies of interest. The fourier transform, however, allows 
one to localize the eigenfrequencies quite accurately (see Table I). 

Near the presumed level crossing, V3 ~ 0.07, the two condensates are fully mixed, see 
Fig. f|. Here, however, the lower frequency, E pa 3.8, is much stronger than the upper one, 
see Fig. |6|. This is clearly seen also at V3 = 0.04 in Fig. [| The amplitudes of the kicked 
condensate, however, depend on the point chosen for the fourier transform. For still stronger 
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couplings t>3 = 0.08, the two condensates are only weakly coupled, see Fig. [|. This may be 
the result of decreasing condensate overlap for large t> 3 ; see Fig. |7|. The amplitudes of the 
oscillations in the two condensates, as shown in Fig. |6|, should not be directly compared; 
only the relative strengths of the oscillational components are indicated. In general, we find 
that level Jj2 is much harder to excite near v 3 m 0.07 than level jjl even if both are easily seen 
in our results. The relative strengths of the two condensate components in one oscillational 
mode varies with position as the corresponding eigen amplitudes. For level Jjl, this is shown 
in Figs. |2H| and f|. We point out that the region where we have found strong mode mixing, 
v\ « t>2 ~ V3, does correspond approximately to the physical situation. This suggests that 
our results may have experimental manifestations in real systems. 
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FIGURES 

FIG. 1. The spectrum of the breathing modes. The lines are the eigenvalues from the Bo- 
goliubov-de Gennes secular equations. The crosses (condensate 2) and the squares (condensate 1) 
is the spectrum calculated with the fourier method. Throughout the calculations we have used 
N\ = N2 = 2100,^1 = 0.02 and V2 = 0.01. All energies are in units of fiV^i^2/2- The levels are 
numbered at V3 = from fjl to $10 starting from below. The fourier method enables us to calculate 
the four lowest modes with a reasonable accuracy, since for higher frequencies the corresponding 
signal decreases rapidly. 

FIG. 2. As an example of our numerical data, this figure shows the spectrum computed by the 
fourier method at the levels ftl and J)2. Here V3 = 0.04 and the two condensates can be seen to 
simultaneously oscillate with the frequency E\ = 3.1 whereas level ft 2 is hardly seen at E% = 4.6. 
The relative amplitude of the oscillations depends on the point where the density is investigated; 
cf. Fig. I 

FIG. 3. The same situation as in Fig.^ with V3 = 0.06. The spectrum shows simultaneous 
oscillations of both condensates at level ftl with E\ = 3.6 and level jj2 with E2 = 4.8. 

FIG. 4. Same situation as in Fig.|2| and Fig.|| with v% = 0.08 where the two condensates are 
decoupled for the two lowest levels. 

FIG. 5. A comparison between the Bogoliubov-de Gennes spectrum and the fourier method 
spectrum calculated with a grid twice as dense as in the previous figures. The general trend with 
a denser grid is that the values approach those from the secular equation, at the expense of a 
fortyfold increase in computing time. 
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FIG. 6. The condensates respond dramatically when uJd approaches a resonance frequency. The 
signals in the left part of the figure are the change in the density as function of time for two different 
driving frequencies. The fourier transform of these signals give the spectra on the right. The time 
is expressed in units of 2/^/^ll^2 and E in l/2h^/^ll^2. Since the interaction strength is here 
V3 = 0.07 the responses of both condensates are significant. Oscillations at frequency J(2 are hardly 
seen here. 

FIG. 7. The more weakly trapped condensate is forming a shell structure separating the con- 
densates with increasing interaction strength V3. Condensate one is concentrated near the center 
of the trap (r ~ 0) and condensate two forms a shell around it. 

FIG. 8. The amplitudes of the eigenmodes for level j)l from the Bogoliubov-de Gennes equations 
shown in a) for and and in b) for and (3^ at V3 = 0.04. 

FIG. 9. The same situation as in Fig. || with V3 = 0.08. The roles of the two condensates are 
changed, which we interpret to derive from passing an avoided crossing. 



17 



TABLES 

TABLE I. A comparison between the Bogoliubov-de Gennes (BdG) solutions and the fourier 
methods. Here Fourieri stands for the crude method whereas Fourier is the method of driving 
the condensate with a frequency ujd- The coalescing oscillations are not shown here. 



i>3 BdG Fourieri Fourier Level (j 



0.002 


3.43 


3.55 


3.50 


1 




4.99 


5.15 


5.10 


2 




6.62 


6.85 


6.60 


3 




9.34 


9.65 


9.60 


4 




9.88 


10.25 


9.80 


5 


0.020 


2.85 


2.95 


3.00 


1 




4.59 


4.80 


4.70 


2 




6.17 


6.35 


6.30 


3 




8.54 


8.90 


9.00 


4 




9.63 


9.90 


10.00 


5 
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